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Interior stagnation point flows of viscoelastic liquids arise in a wide variety of appli- 
cations including extensional viscometry, polymer processing and microfluidics. Experi- 
mentally, these flows have long been known to exhibit instabilities, but the mechanisms 
underlying them have not previously been elucidated. We computationally demonstrate 
the existence of a supercritical oscillatory instability of low-Reynolds number viscoelastic 
flow in a two-dimensional cross-slot geometry. The fluctuations are closely associated with 
the "birefringent strand" of highly stretched polymer chains associated with the outflow 
from the stagnation point at high Weissenberg number. Additionally, we describe the 
mechanism of instability, which arises from the coupling of flow with extensional stresses 
and their steep gradients in the stagnation point region. 



1. Introduction 

While Newtonian flows become unstable only at high Reynolds number Re, when the 

inertial terms in momentum balance dominate, flows of viscoelastic fluids such as polymer 

solutions and melts are known to have interesting instabilities and nonlinear dynamical 
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behaviors even at extremely low Re. These "purely e lastic" instab i lities arise in rh eometry 



of complex fluids as well as in many applications ([Larson 



f992; 



Shaqfeh 



1996). Recent 



studies of viscoelastic flows in microfluidic d evices broaden the scop e of these nonlin- 



2005). The small length 



ear dynamical problems of viscoelastic flows (jSauires fc Quakd ll 
scales in microfluidic devices enable large shear rates, and thus high Wi (Weissenberg 
number, Wi = Xy, where A is a characteristic time scale of the fluid and 7 is a charac- 
teristic shear rate of the flow), at very low Re. Instabilities are not always undesirable, 
especially when the accompanying flow modification is controllable and can thus be uti- 
lized in the design and operation of microfluidic devices. Specifically, instabilities have 
been found and flow-controlling logic elements have been desi gned in a series of microflu - 



idic geometries, e. g. flow rectifier with a nisotropic resistance (|Groisman fc Q uake 



flip-fl op memory (jGroisman et al 



2003) and nonlinear flow resistance 



2004), 



Groisman et al 



20031) . Another prospective appl ication of these instabilities is to enhancement of mixing 



at lab-on-a-chip length scales (jGroisman fc Steinberg 



2001), where turbulent mixing is 



absent due to small length scales and an alternative is needed. 

The best understood of these instabilities are those that occur in viscometric flows with 



curve d streamlines: e.g . flows in Taylor-Couett e (jMuller et 



1994) 



2000: 



Magda fc Larson 



, cone-and-plat e ( Magda fc Larsonll988l ) and parallel-plates (jGroisman fc Steinberg 



1989 ). Tavlor-Dean Ooo fc 



Shaqfeh 



1988) flow geometries. In these geometries, the primary source of 



instability is the coupling of normal stresses with streamline curvature (i.e. the pr esence of 



"hoop stresses"), leading to radial compressive forces that can drive instabilities flSh aq fch 



199 



Muller et al 



198£ 



Pakdel fc McKinlev 



Joo fc Shaqfeh 



199' 



Graham 



coelastic free surface flows (jSpiegelberg fc McKinlev 



1994; 



Magda fc Larson 



1988 



Larson et al 



1990; 



1998). Similar mechanisms drive instabilities in vis- 



1996 



Graham 



2003). 



Attention in this paper focuses on a different class of flows, whose instabilities are not 



Viscoelastic cross-slot flow instability 
well-understood - stagna tion point fl ows, like those gener ated with oppo sed-iet (IChow et al. 



1988; 



Muller et al. 



1988 ), cross-slot ( Arratia et al 



and four-roll mill (|Ng fc Leal 



1993 



Broadbent et al 



2006), two-roll mill ( Ng fc Lea 



1993) 



19781) devices. Figure Q] shows a 



schematic of a cross-slot geometry. A characteristic phenomenon in these stagnation 
point flows is the formation of a narrow region of fluid with high polymer stress ex- 
tending downstream from the stagnation point. This region can be obser ved in optical 



exper iments as a bright bire fringent "strand" with the rest of flu id dark (jHarlen et 



1990). Keller and coworkers ( Chow et 



1988 



Muller et al. 



1988) reported instabilities 



in stagnation point flows of semi-dilute polymer solutions generated by an axisymmetric 
opposed-jet device. Specifically, for a fixed polymer species and concentration, upon a 
critical extension rate (or critical Wi) polymer chains become stretched by flow near the 
stagnation point and a sharp uniform birefringent stand forms. The width of this birefrin- 
gent strand increases with increasing Wi until a stability limit is reached, beyond which 
the birefringent strand becomes destabilized and changes in its morphology are observed. 
At higher Wi, the flow pattern and birefringent strand become time-dependent. Recent 
tracer and particle-tracking experimen ts of stagnation poin t flow in a micro-fabricated 



20061 ) show instabilities of dilute 



cross-slot geometry by Arratia et al. (jArratia et 
polymer solution at low Re (< 10~ 2 ). In their experiments fluid from one of the two 
incoming channels is dyed and a sharp and flat interface between dyed and undyed fluids 
is observed at low Wi. Upon an onset value of Wi, this flow pattern loses its stability: 
spatial symmetry is broken but the flow remains steady. The interface becomes distorted 
in such a way that more than half of the dyed fluid goes to one of the outgoing channels 
while more undyed fluid travels through the other. At even higher Wi the flow becomes 
time-dependent and the direction of asymmetry flips between two outgoing channels with 
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time. Particle-tracking images in the time-dependent flow pattern indicate the existence 
of vortical structures around stagnation point. 

Another class of stagnation point flows is associated with liquid-solid or liquid-gas 
interfaces, such as flows passing submerged solid obs tacles, around moving bubbles or 



19931 ) reported a 



toward a free surface. For example, McKinley et al. (jMcKinlev et al. 
three-dimensional steady cellular disturbances in the wake of a cylinder submerged in 
a viscoelastic fluid. Around a falling sphere in viscoelastic fluids, fore-and-aft symmetry 
of velocity field is broken and the velocity perturbation in the wake can be away from 
the sphe re, toward the s phere or a combination of the two depend ing on the polymer 



solution ( Hassager 



1979; 



Remmelgas et al. (Rc mmelgas et 



Bisgaard & Hassager 



1982; 



Bisgaard 



1983). 



19991 ) computationally studied the stagnation 



point flow in a cross-slot geometry with two different FENE (finite extensible nonlinear 
elastic) dumbbell models. Using the two models, they studied the effects of configuration- 
dependent friction coefficient on polymer relaxation and the shape of the birefringent 
strand. Their simulation approach was restricted to rela tively low W i (<~ O(l)) with 



20021 ) conducted sim- 



symmetry imposed on centerlines of all channels. Harlen (jHarlenl 
ulations of a sedimenting sphere in a viscoelastic fluid to explore the wake behaviors. 
He explains the experimental observations of both negative (velocity perturbation away 
from the sphere) and extended (velocity perturbation toward the sphere) wakes in terms 
of combined effects of the stretched polymer in the birefringent strand following the 
stagnation point behind the sphere and the recoil outside of the strand. Neither of these 
analyses directly addressed instabilities of these flows. 

Various approximate approaches have been taken in the past to obt ain an understand- 



ing of these instabilities observed in experiments. Harris and Rallison ([Harris &; Rallison 



1993, 



19941 ) investigated the instabilities of the birefringent strand behind a free isolated 
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stagnating point through a simplified approach, in which polymer molecules are modeled 

as linear-locked dumbbells, which are fully stretched within a thin strand lying along the 
centerline. Polymer molecules contribute a normal stress proportional to the extension 
rate only when they are fully stretched (i.e. in the strand), otherwise the flow is treated 
as Newtonian. The lubrication approximation is applied for the Newtonian region and 
the effects of birefringent strand are coupled into the problem through point forces along 
the strand. Two instabilities are reported. At low Wi (« 1.2 — 1.7), a varicose disturbance 
is linearly unstable, in which the width of birefringent strand oscillates without breaking 
the symmetry of the flow pattern. At higher Wi another instability is observed in which 
symmetry with respect to the extension axis breaks and the birefringent strand becomes 
sinuous in shape and oscillatory with time, with zero displacement at the stagnation point 
and increasing magnitude of displacement downstream from it. Symmetry with respect 
to the inflow axis is always imposed. The mechanism of these instabilities is explained: 
perturbations in the shape or position of the birefringent strand affect the stretching of 
incoming polymer molecules such that they enhance the perturbation after they become 
fully stretched and merge into the strand. This mechanism is close to the one we are 
about to present later in this paper with regard to the importance of flow kinematics 
and the extensional stress. However, in their linear stability analysis the spatial depen- 
dence of the birefringent strand in the outflow direction is neglecte d, which is important 



19971 ) studied steady state 



according to our simulations. Oztekin et al. (jOztekin et 
similarity solution for planar stagnation point flow at a solid wall predicting that this 
flow is linearly unstable to local three-dimensional disturbances. Their results indicate 
that traveling wave disturbances that have periodic structure in the neutral direction 
could lead to instabilities of steady state solutions above certain critical Wi. 

In this paper, we present numerical simulation results of viscoelastic stagnation point 
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flow in a two-dimensional cross-slot geometry. With increasing Wi, we observe the forma- 
tion and elongation of the birefringent strand across the stagnation point. At high Wi, we 
find the occurrence of an oscillatory instability. These results re semble the experi mental 



observations of oscillatory birefringent width by Miiller et al. ( Miiller et al 



1988) and 



1994J). By analyzing 



the varicose instability predicted by Harris et al. (jHarris fc Rallisonl l 
the perturbations in both velocity and stress fields, a novel instability mechanism based 
on normal stress effects and flow kinematics is identified. 



2. Formulation and Methods 

We consider a fourfold symmetric planar cross-slot geometry, as shown in Figure [TJ 
Flow enters from top and bottom and leaves from left and right. For laminar Newtonian 
flow, two incoming streams meet at the intersection of the cross and each of them splits 
evenly and goes into both outgoing channels, generating a stagnation point at the origin 
near which an extensional flow exists. We use round corners at the intersections of channel 
walls in order to avoid enormous stress gradients at the corners, which cause numerical 
difficulties. 

The momentum and mass balances are: 

Re(— + u- Vu \ = -Vp + PV 2 u + (1 - j3) — (V • r p ) , (2.1) 

V-m = 0. (2.2) 

Parameters in Equations (|2.ip and (|2.2[) are defined as: Re = pUl/ (r] s + rj p ), Wi = 2XU/1 
and (3 = ry s / (?7 S + ri p ), where p is the fluid density and for dilute polymer solution we 
assume it to be the same as the solvent density; rj s is the solvent viscosity and rj p is 
the polymer contribution to the shear viscosity at zero shear rate and U and I are 
characteristic velocity and length scales of the flow. Here I is chosen to be the half 



Viscoelastic cross-slot flow instability 



7 



21 



10/ 



^^^^^ 



0.5/ 



Figure 1. Schematic of the cross-slot flow geometry, 
channel width and the definition of U is based on the pressure drop applied between the 
entrances and exits of the channel. Specifically, U is defined to be the centerline velocity 
of a Newtonian plane Poiseuille flow under the same pressure drop in a straight channel 
with length 201, which is comparable to the lengths of streamlines in the present geometry. 
According to this definition, the nondimensional pressure drop in our simulation is fixed 
at 40 and the centerline Newtonian velocity in cross-slot geometry is typically slightly 
lower than 1 since the extensional flow near the stagnation point has a higher resistance 
than that in a straight channel. The polymer contribution to t he stress tensor is denoted 



t p and is calculated with the FENE-P constitutive equations ([Bird et 



1987): 



a Wi (da _ _ . _ , T 

- . , -I — — h u ■ Va — a ■ Vit — (a • Vw 

. tr(oQ 2 V dt 



6 + 5 



a 



b \ i_ Msd V 6 + 2 

\ b 



6 + 2 
2 



(2.3) 
(2.4) 



In Equations (|2.3[) and (|2.4[) . polymer chains are modeled as FENE dumbbells (two beads 
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connected by a finitely-extensible- nonlinear-elastic spring). Here a = (QQ) is the con- 
formation tensor of the dumbbells where Q is the end-to-end vector of the dumbbells and 
{•) represents an ensemble average. The parameter b determines the maximum extension 
of dumbbells, i.e. the upper limit of tr(cc). 

At the entrances and exits of the flow geometry, normal flow boundary conditions are 
applied, i.e. t ■ u = where t is the unit vector tangential to the boundary. Pressure is 
set to be 40 at entrances and at exits. No-slip boundary conditions are applied at all 
other boundaries. Boundary conditions for stress are only needed at the entrances, where 
the profile of a is set to be the same as that for a fully developed pressure-driven flow 
in a straight channel with the same Wi. Other fixed parameters in our simulations are: 
Re = 0.1, (3 = 0.95 and b = 1000, which means we focus on dilute solutions of long-chain 
polymers at low Reynolds number. 



The discrete elastic stress splitting (DEVSS) formulation (jBaaiiens el al. 



1997 



Baaiiens 



19981 ) is applied in our simulation: i.e. a new variable A is introduced as the rate of strain 



and a new equation is added into the equation system: 



A = Vu + V/. (2.5) 



A numerical stabilization term 7V • ^Vit + Vu T — Aj is added to the right-hand-side 
of the momentum balance (Equation (|2.ip ) and it is worthwhile to point out that this 
term is only nontrivial in the discretized formulation and does not change the physical 
problem. In this term, 7 is an adjustable parameter and 7 = 1.0 is used in our sim- 
ulations. The velocity field u is interpolated with quadratic elements while pressure p, 
polymer conformation tensor a. and rate of strain A are interpolated with linear ele- 



ments. Consistent with Baaijens's conclusion (jBaaiiens 



19981 ). DEVSS greatly increases 



the upper limit of Wi achievable in our simulations. Quadrilateral elements are used for 
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all variables. Our experience shows that quadrilateral elements have great advantages 

over triangular ones, yielding much better spatial smoothness in the stress field at com- 
parable degrees of freedom to be solved. Another merit of quadrilateral elements is the 
capability of manual control over mesh grids. This is extremely important when certain 
restrictions, such as symmetry, are required. In our simulation, finer meshes are used 
within and around the intersection region of the geometry and the mesh is required to 
be symmetric with respect to both axes. Within a horizontal band (—0.2 < y < 0.2) 
across the stagnation point, very fine meshes are generated to capture the sharp stress 
gradient along the birefringent str and. The streamline upwind/Pctrov-Galerkin method 



(SUPG) ferooks fc Hughes 



1982) is applied in Equation (|2.3p by replacing the usual 
Galerkin weighting function w with w + Shu ■ Vw/||it||, where h is the geometric average 
of the local mesh length scales and S is an adjustable parameter, set to 8 = 0.3 in our 
simulations. This formulation is implemented using the commercially available COMSOL 
Multiphysics software. 

3. Results and Discussions 

3.1. Steady States 

Steady state solutions are found for all Wi investigated (0.2 < Wi < 100) in our study. 
For Wi ^ 60 steady states are found by time integration and for those with larger Wi 
Newton iteration (parameter continuation) is used because of possible loss of stability, 
as we describe below. At low Wi the velocity field is virtually unaffected by the polymer 
molecules. Velocity contours at Wi — 0.2 are plotted in Figure 2(a)[ for clarity only part 



of the channel is shown. A stagnation point is found at the center of the domain ((0, 0)) In 
both incoming and outgoing channels, the flow is almost the same as pressure driven flow 
in a straight channel. No distinct difference can be observed for the incoming and outgoing 
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directions in velocity field. Figure 2(b) shows contours of extension rate at Wi = 0.2, 
in which a region dominated by extensional flow is found near stagnation point. High 
extension rate is also found near the corners due to the no-slip walls. The magnitude of 
polymer stretching can be measured by the trace of its conformation tensor tr(a), and is 



plotted in Figure 2(c) At low Wi, the extent to which polymers are deformed is barely 
noticeable, but it can be clearly seen that polymers are primarily stretched in either the 
extensional flow near the stagnation point and corners or the shear flows near the walls. 
At high Wi (Wi = 50, Figure [3]), the situation is very different. Polymers are strongly 
stretched by the extensional flow near the stagnation point and this stretching effect by 
extensional flow overwhelms that of the shear flow. A distinct band of highly stretched 



polymers (the birefringent strand) forms (Figure 3(c) ). Since the polymer relaxation time 
in this case is larger than the flow convection time from stagnation point to the exits, this 
birefringent strand extends the whole length of the simulation domain. The resulting high 



polymer stress significantly affects the velocity field (Figure 3(a) ). Regions with reduced 
velocity extend much farther away in the downstream directions of the stagnation point 
than in the low Wi case, especially along the x-axis, where high polymer stress dominates. 
Correspondingly, a reduction in the extension rate near the stagnation point is observed, 
most noticeably along the birefringent strand (Figure |3(b)| . 

Figures [4] and [5] show profiles at various values of Wi of tr(a) along the outflow (x-axis) 
and inflow (y-axis) directions of this stagnation point (note the difference in scales in the 
two plots). For increasing Wi the length of the region with highly stretched polymer keeps 
increasing due to the increased relative relaxation time (Figured]). In high Wi cases (Wi = 
30 and Wi = 100), polymers are not fully relaxed even when they reach the exit of the 
simulation domain. The cross-sectional view of tr(a) profile along the y-axis (Figure [5]) 
shows interesting non-monotonic behaviors. Although the height of the profile (tr(a) max ) 
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(a) H| 




-0.6 -0.4 -0.2 0.2 0.4 0.6 0.£ 
(b) dux/dx 




2.99 3.01 3.03 3.05 3.07 3.09 3.11 3.13 
(c) tr(a) 

Figure 2. Contour plots of steady state solution: Wi = 0.2 (only the central part of the flow 

domain is shown). 

keeps increasing upon increasing Wi, the width of Wi = 100 case is smaller than that 
of Wi — 30, resulting in a steeper transition section between low and high stretching 
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0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 
(a) H| 




-0.6 -0.4 -0.2 0.2 0.4 0.6 0. 
(b) e = dux/dx 




100 200 300 400 500 600 700 800 900 
(c) tr(a) 

Figure 3. Contour plots of steady state solution: Wi = 50 (only the central part of the flow 

domain is shown). 

regions. If we arbitrarily define tr(a) > 300 as the observable birefringence region, the 
width W and the length L of the birefringent strand (measured on the inflow and outflow 
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Figure 4. Profile of tr(a) along y = 0. 

axes, respectively) can be plotted as functions of Wi, as in Figure [6] (values of L for 
Wi > 30 are not shown since they exceed the length of the simulation domain) . A clear 
non-monotonic trend is observed in the plot of birefringence width, where W increases 
sharply at relatively low Wi and peaks around Wi — 40. After that W decreases mildly 
but consistently with further higher Wi. This non-monotonic tr end is consistent w ith 



1988) 



experimental observations of birefringence in opposed-jet devices ([Miiller et al. 

Similarly, a non-monotonicity is also found in the change of velocity field with Wi. 
Shown in Figure [7] is the value of extension rate, averaged within a box around the 
stagnation point (—0.1 < x < 0.1, —0.1 < y < 0.1), as a function of Wi. As Wi increases, 
the extension rate decreases at low Wi but increases at high Wi, with a minimum found 
around Wi — 40. Besides, most of experimental results are presented in terms of Deborah 
number (De), defined as the product of the polymer relaxation time and an estimate of the 



extension rate near the stagnation point. Noticing that the average (nondimensionlized) 
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Figure 5. Profile of tr(a) along x = in the region very near the stagnation point. 
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Figure 6. Birefringence strand width W; Inset: Birefringence strand length L (tr(a) > 300 is 
considered as observable birefringence region). 

extension rate changes within a very narrow range (around 0.55 ~ 0.6), a conversion 
De = Q.iWi can be adopted for comparison of our results with experimental ones. 
Some understanding of this non-monotonicity can be gained by looking at Figure 0J 
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Figure 7. Average extension rate (du x /dx) avg (averages taken in the domain 
-0.1 < x < 0.1, -0.1 < y < 0.1). 

Here it can be seen that for Wi < 30, the birefringent strand is not yet "fully developed" 
in the sense that the polymer stretching is not yet saturating near full extension. Thus 
the evolution of the velocity field in this regime of Wi reflects the significant changes 
that occur in the stress field in this regime. At higher Wi, however, the polymer stress 
field in the strand is saturating, and thus not changing significantly. Furthermore, at 
these high Wcisscnbcrg numbers, the relaxation of stress downstream of the stagnation 
point diminishes, decreasing the gradient dr xx /dx and thus decreasing the effect of 
viscoelasticity on the flow near the stagnation point. 

3.2. Periodic Orbits 

We turn now to the stability of the steady states that have just been described. Rather 
than attempting to compute the eigenspectra of the linearization of the problem, an 
exceedingly demanding task, we examine stability by direct time integration of perturbed 
steady states. The perturbations take the form of slightly asymmetric pressure profiles at 
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0.278-, 




0.262 - 
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0.156 0.158 0.160 0.162 0.164 

w 

Figure 8. Two dimensional projection of the dynamic trajectory from the steady state to the 
periodic orbit at Wi = 66: ||u|| = \Ju% + u\ at (0.5,0) v.s. W. 

the two entrances (0.1% maximum deviation from the steady state value) that are applied 
for one time unit, then released. As an example, Figure [8] shows a two dimensional 
projection of the trajectory of the system evolution over time at Wi — 66. Here the 
velocity magnitude at a point near the stagnation point ((0.5,0)) is plotted against the 
bircfringent strand width W measured on the inflow axis. The system starts at the steady 
state with W = 0.1593 and ||u|| (0.5,0) = 0.2687 and spirals outward with time after the 
perturbation. Eventually the trajectory merges into a cycle (the outer dark cycle in 
the Figured]). This clearly identifies the existence of a stable periodic orbit. Note the 
anticorrelation between ||u|| and W, i.e. when the flow speeds up near the stagnation 
point, the strand thins and vice versa. 

Figure [5] shows the root-mean-square deviations over one period of W from its steady 
state values , normalized by the corresponding steady state values W s . s ., as a function of 
Wi for all the cases where we found periodic orbits. Time integrations for Wi > 74 did 
not converge due to the enormous stress gradient around the corners of the no-slip walls 
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Figure 9. Left: Root mean square deviations of the birefringent strand width W at periodic 
orbits, normalized by steady state values; Right: Oscillation periods. 

and the consequent numerical oscillations downstream. Data points for W rms computed 
from our simulations are fitted with a function of the form a(Wi — b) c , with c fixed 
at 1/2. Very good agreement is found for our sim ulation data with the 1/2 powe r law, 



1983). The 



characteristic of a supercritical Hopf bifurcation (jGuckenheimer fc Holmesl 
critical Weissenberg number Witriticea is identified to be 64.99 by this fitting. Also shown 
in Figure O are periods of oscillations, where a slight decrease with increasing Wi is 
found. This is interesting since it indicates that some time scale other than the polymer 
relaxation time sets the period of oscillations. 

Time-dependent fluctuatio ns of birefringence width are also reported in the experi- 



ments done by Miiller et al. (jMuller et al 



1988) in opposed-jet devices. In their optical 



experiments with semi-dilute aPS solutions, the width of the birefringent strand oscil- 
lates rapidly between two values in a certain range of extension rate. The critical value of 
extension rate for the instability in their study is close to the one where the birefringent 
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strand width W is highest, while in our simulations Wicritic&l (— 65) is somewhat larger 
than the one ( Wi w 40) that gives the largest W. 

3.3. Instability Mechanism 

We turn now to the spatiotemporal structure of the instability and its underlying physical 
mechanism. We will denote the deviations in velocity, pressure and stress with primes, 
while steady state values will be denoted with a superscript "s" : 

u = u s + u', (3.1) 
P = P s +p', (3.2) 
a = a s + a'. (3.3) 

Figures [TOl [TT1 and [T2l illustrate u' x , u' y and ot' xx , respectively, at intervals of 1/8 period, 
corresponding to the periodic orbit at a Weissenberg number close to the bifurcation 
point ( Wi = 66). Time starts from an arbitrarily chosen snapshot on the periodic orbit 
and only a quarter of the region near the stagnation point is shown, behavior in the rest 
of the domain can be inferred from the reflection symmetry across the axes. 



At the beginning of the cycle (Figure 10(a)), u' x is positive in the region very close 



to the stagnation point while it is negative in most of the downstream region. As time 
goes on, this positive deviation near the stagnation point grows into a "jet" , a region of 
liquid moving downstream away from the stagnation point faster than the steady state 
velocity, as shown in Figures [lO(b)[ [l0(c") and |10(d")| Correspondingly, by continuity, the 



inflow toward the stagnation point is also faster as shown in Figures 11(a) -11(d) Note 



that very near the stagnation point deviations from steady state remain small. At the 



beginning of the second half of the cycle (Figure 10(e)), the jet extends further down 



stream and grows to the full width of the channel. Meanwhile, in the region closer to the 
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stagnation point, velocity deviations drop (Figures 10(e) 11(e) ) and start to change signs 
(Figures |10(f)[ [Tl (f ) ). Consequently, the growth of the jet is interrupted and a "wake", 
a region of fluid moving slower than the steady state velocity, emerges downstream (Fig- 
ures |10(f)] - |10(h)| and 11(f) - 11(h) ). Similarly, as the wake grow larger, velocity deviations 



near the stagnation point change signs and a new cycle starts (Figures 10(a) and 11(a)). 

The velocity deviations are closely related with those of the stress field (Figure [T2|). 
Generally speaking, "jets" are accompanied by negative a' xx and thus thinning of the 
birefringent strand and "wakes" are associated with the birefringent thickening. The 
largest deviations are found at the edges of the birefringent strand where da s xx /dy is 
largest. Note that deviations in the stress field are always small along the center line of 
the birefringent strand because there polymer molecules are almost fully stretched and 
the huge spring force is sufficient to resist any perturbations. 

One may notice the small spatial oscillations in the stress field deviations, characterized 
by alternating high and low stress stripes, along the outflow direction. These oscillations, 
apparently unphysical and centered around zero, also exist along the birefringence strand 



in steady state solutions, though they are not easy to see from the contour in Figure 3(c) 
as they are overwhe lmed by high tr(ct) in the birefringent strand. Unfortunately, as 



proven by Renardy (jRenardv 



2006), spatial non-smoothness is inevitable in numerical 
simulations of viscoelastic extensional flow upon certain Wi due to the singularities in 
stress gradients. These singularities could not be fully resolved by any finite mesh size 
and this problem would always show up in numerical solutions of high Wi viscoelastic 
stagnation point flows. However, we do not expect these oscillations to qualitatively 
affect our observations for a couple of reasons. First, non-smoothness has been observed 
in our simulation at Wi values much lower than the critical Wi of this instability. Second, 



observable non-smoothness is always found some distance away from the stagnation point 



20 L. Xi AND M. D. Graham 

in the downstream direction while the instability is dominated by the physics in the close 
vicinity of the stagnation point and since FENE-P is a convective equation we do not 
expect anything occurring downstream to affect upstream dynamics. Last, and most 
importantly simulations with different meshes display different mesh size dependent 
stripes, while the nature of the instability remains virtually unchanged. 

Insight into the mechanism of this instability can be gained by examining the linearized 
equation for a' xx : 

da' xx _ 2 a' xx 2 < x tr(a') 



dt Wll- t Afl Wi b ^ _ tr(a') ^ 



2 



, s da 'xx s da 'xx , 9a s xx , da s xx (3.4) 
x dx y dy x dx y 8y 

+ za xx Qx + za xy Qy + za xx ^ + za xy ^ . 

In the following analysis, terms on the right-hand-side (RHS) of Equation l3.4l are named 
"RHS*", where "*" is determined by the order of appearance on the RHS. Terms and 
their physical meanings are summarized in Table [TJ To understand the mechanism of the 
instability, magnitudes of these terms at the point (0, —0.05) are plotted as a function of 
time during roughly a period in the bottom view of Figure[T51 Terms RHS3, RHS5, RHS8 
and RHSIO are zero by symmetry and not plotted. This position is right at the edge of 
the birefringent strand and as shown in Figure [T2l it is also where significant deviations 
in the stress field are observed. Time-dependent oscillations at other places, including off 
the symmetry axis x = 0, have also been checked and nothing that could qualitatively 
affect our analysis was seen. Correspondingly, deviations in polymer conformation, inflow 
velocity and extension rate, normalized by steady state values, are plotted in the top view 
of Figure [131 

Consistent with our earlier observations, deviations in the velocity field [u' y and du' x /dx) 
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(h) £ = 11.74 



Figure 10. Perturbation of x-component of velocity, u' x with respect to steady states at periodic 
orbits: Wi = 66. The region shown is < x < 1.5, —1.1 < y < 0, stagnation point is at the 
top-left corner. 

and deviations in stress field (a' xx ) are opposite in signs for most of the time within the 
period. Among the terms plotted, RHS4, RHS6, RHS7 and RHS9 are much larger than 
the relaxation terms, RHS1 and RHS2, and dominate the dynamics. (Relaxation terms 
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(a) t = (b) t = 1.68 



C.b .... !. ... I- = C.b = 




(c) t = 3.35 (d) t = 5.03 



1.5 C C.b .... !. ... 




(g) t = 10.06 (h) t = 11.74 



Figure 1 1 . Perturbation of y-component of velocity, u' y with respect to steady states at periodic 
orbits: Wi = 66. The region shown is < x < 1.5, —1.1 < y < 0, stagnation point is at the 
top-left corner. 

are large at the very inner regions of the birefringcnt strand and that is why oscillations 
in the stress field there are barely noticeable.) Moreover, RHS4, RHS6 and RHS9 are 
mostly in phase with a' xx and thus tend to enhance the deviations while RHS7 is out of 



Viscoelastic cross-slot flow instability 23 




(g) t = 10.06 (h) t = 11.74 



Figure 12. Perturbation of xx-component of polymer conformation tensor, a' xx with respect 
to steady states at periodic orbits: Wi = 66. The region shown is < x < 1.5, —1.1 < y < 0, 
stagnation point is at the top-left corner. The edge of the steady state birefringent strand is the 
line y « —0.05. 

phase with a' xx and hence damps the deviations. It is the joint effect of these compet- 



ing destabilizing and stabilizing forces that gives the oscillatory behavior of the system. 
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Term 


Formula 


Physical Significance 


RHS1 


2 

Wi , tr(a») 
1 b 


Relaxation. 


RHS2 








j Relaxation. 


RHS3 




Convection of conformation deviations by the 
steady state x-velocity. 


RHS4 


„,s S a xx 


Convection of conformation deviations by the 
steady state y-velocity. 


RHS5 


/ a «xx 


Convection of the steady state conformation by 
x-velocity deviations. 


RHS6 


/ 9a| x 
U V dy 


Convection of the steady state conformation by 
y-velocity deviations. 


RHS7 


a "' c dx 


Stretching caused by deviations in the extension 
rate. 


RHS8 


2a" ^ 


Stretching caused by deviations in the shear rate. 


RHS9 




Stretching caused by deviations in the extensional 
stress. 


RHS10 2a4 s ^ 


Stretching caused by deviations in the shear 
stress. 



Table 1. Terms on the right-hand side of Equation (|3.4 



Finally, notice that among the three destabilizing terms, RHS6 is the one that leads the 
phase and thus guides the instability. 

Based on these observations from Figure [131 a mechanism for the instability can be 
proposed, which is illustrated schematically in Figure [TH At the beginning of the cycle 
(t = 0), u' y is slightly above zero, indicating that the inflow speed is faster than that 
in the steady state. As a consequence, RHS6 becomes negative first, followed by RHS4 
and RHS9. In particular, a faster incoming convective flow brings unstretched polymer 
molecules toward the stagnation point (corresponding to RHS6), as depicted in Fig- 



ure 14(a) These polymer chains have less time to get stretched and when they reach the 
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Figure 13. Time-dependent oscillations at (0, —0.05). Top view: perturbations of variables nor- 
malized by steady-state quantities; Bottom view: magnitudes of terms on RHS of Equation (|3.4[) . 



edges of the birefringent strand (e.g. dumbbell B), they are less stretched compared with 
the steady state. As a result, fluid around dumbbell B has lower stress than at the steady 
state, corresponding to a thinning of the birefringent strand. Meanwhile, since dumbbell 
B contains smaller spring forces than its downstream neighbors A and A', the net forces 
(shaded arrows) exerted by polymer on the fluid point outward, generating jets down- 
stream from the stagnation point. (In other words, when the stress at the center is lower, 
the net stress divergence points outward, which increases momentum in the downstream 
directions.) By continuity, more fluid has to be drawn toward the stagnation point and 
the initial deviation in u' y is then enhanced. However, as the flow speeds up in the vicinity 
of the stagnation point, the extension rate also starts to increase. This effect (correspond- 
ing to RHS7) tends to stretch polymer molecules more and stabilize the deviations, as 
shown in Figure [T31 Eventually this effect will be able to overcome that of RHS6 as well 
as RHS4 and RHS9 and the stress near the stagnation point starts increase after it passes 
the minimum at around t = 3.5, which causes a re-thickening of the birefringent strand 
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Centerline 




(a) Thinning process of the birefringent strand. 



Centerline 




(b) Re-thickening process of the birefringent strand. 



Figure 14. Schematic of instability mechanism (view of the lower half geometry) (Gray solid 
arrows: spring forces of the dumbbells; Large shaded arrows: net forces exerted by polymer 
molecules (dumbbells) on the fluid). 



as illustrated in Figure |14(b)| By a similar argument as that above, dumbbell C has 



higher spring forces than B and B', the dumbbells which were passing near the center 
when stress was at minimum, and the net polymer forces point inward, which starts to 
suppress the jets. Inflow velocity decreases as the birefringent strand thickens and this 
gives incoming polymer molecules more time to be stretched and further thickens the 
birefringent strand. Eventually a xx will come back to the steady state value at around 
t = 7.2, however, since all the deviations are not synchronized, a negative deviation is 
found in u v and an identical analysis with opposite signs can be made for the second half 
of the cycle. 

Within this mechanism, a sharp edge of the birefringent strand, i.e. large magnitude 
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of da xx /dy (~ 0(1O 4 ) in our simulations), is required so that a small u' y can give a 

sufficiently large RHS6 to drive the instability. This is made possible by the kinemat- 
ics of the flow near the stagnation point, where the incoming polymer molecules are 
strongly stretched within a short distance. Another similar effect is that stress deriva- 
tives are stretched in the outgoing direction and thus greatly weakened as fluid moves 
downstream; therefore the instability is dominated by physics in the vicinity of the stag- 
nation point. In the earlier mechanism for th e so-called "varicose instability", given by 



Harris and Rallison (Harris &: Rallison 



1994 ) . the importance of extensional stress and 



flow kinematics, especially the role of convection of incoming molecules, was also recog- 
nized. However, the picture described in their work is not the same as ours due to the 
simplifications in their model. Their linear stability analysis ignores the x-dependence 
of the birefringent width while in our simulations, x-dependence of the stress field is 
closely related to the changes in velocity field. Besides, their analysis does not identify a 
restoring force for the deviations and the oscillatory behavior could not be explained. 



4. Conclusions 

Using a DEVSS/SUPG formulation of the finite clement method, we are able to simu- 
late viscoelastic stagnation point flow and obtain steady state and time-dependent solu- 
tions at high Wi. For Wi 3> 1, a clear birefringent strand is observed. The width of this 
birefringent strand increases with increasing Wi until Wi ~ 40 after which it declines 
gradually. This also results in a non-monotonic trend in the modification of the velocity 
field. 

At around Wi = 65 the steady state solution loses stability and a periodic orbit be- 
comes the attractor in phase space. Flow motion of the periodic orbit is characterized by 
time-dependent fluctuations, specifically, alternating positive (jet) and negative (wake) 
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deviations from the steady state velocity in the regions downstream of the stagnation 
point. A mechanism is proposed which, taking account of the interaction between ve- 
locity and stress fields, is able to explain the whole process of the oscillatory instability. 
Extensional stresses and their gradients as well as the flow kinetics near the stagnation 
points is identified as important factors in the mechanism. This mechanism is different 
from that of the "hoop stress" instabilities, which occur in viscometric flows with curved 
streamlines, and we expect that this mechanism could be extended and explain various 
instabilities occurring in viscoelastic flows with stagnation points. 
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